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Abstract 


The coupling of molecular dynamics (MD) simulations with finite 
element methods (FEM) yields computationally efficient models that link 
fundamental material processes at the atomistic level with continuum field 
responses at higher length scales. The theoretical challenge involves 
developing a seamless connection along an interface between two 
inherently different simulation frameworks. Various specialized methods 
have been developed to solve particular classes of problems. Many of 
these methods link the kinematics of individual MD atoms with FEM 
nodes at their common interface, necessarily requiring that the finite 
element mesh be refined to atomic resolution. Some of these coupling 
approaches also require simulations to be carried out at 0 K and restrict 
modeling to two-dimensional material domains due to difficulties in 
simulating full three-dimensional material processes. In the present work, 
a new approach to MD-FEM coupling is developed based on a 
restatement of the standard boundary value problem used to define a 
coupled domain. The method replaces a direct linkage of individual MD 
atoms and finite element (FE) nodes with a statistical averaging of 
atomistic displacements in local atomic volumes associated with each FE 
node in an interface region. The FEM and MD computational systems are 
effectively independent and communicate only through an iterative update 
of their boundary conditions. With the use of statistical averages of the 
atomistic quantities to couple the two computational schemes, the 
developed approach is referred to as an embedded statistical coupling 
method (ESCM). ESCM provides an enhanced coupling methodology that 
is inherently applicable to three-dimensional domains, avoids 
discretization of the continuum model to atomic scale resolution, and 
permits finite temperature states to be applied. 
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1. Introduction 


The emerging field of nanomechanics is providing a new focus in the study of the 
mechanics of materials, particularly that of simulating fundamental atomic mechanisms 
involved in the initiation and evolution of damage. These simulations are commonly 
based on either quantum mechanics (ab-initio, tight-binding (TB) or density-functional 
theory (DFT)) methods, on classical molecular dynamics (MD), or molecular statics 
(MS) methods. These predictions of material behavior at nanometer length scales promise 
the development of physics-based ’bottom-up’ multiscale analyses that can aid in 
understanding the evolution of failure mechanisms across length scales. However, 
modeling atomistic processes quickly becomes computationally intractable as the system 
size increases. With current computer technology, the computational demands of 
modeling suitable domain sizes (on the order of hundreds of atoms for quantum 
mechanics-based methods, and potentially billions of atoms for classical mechanics- 
based methods) and integrating the governing equations of state over sufficiently long 
time intervals, quickly reaches an upper bound for practical analyses. In contrast, 
continuum mechanics methods such as the finite element method (FEM) provide an 
economical numerical representation of material behavior at length scales in which 
continuum assumptions apply. However, all constitutive relationships, kinematics, etc, 
must be assumed a priori. 

The concept of bridging length scales by concurrently coupling atomistic and 
continuum computational paradigms is particularly attractive as a highly efficient means 
of reducing the computational cost of simulations in cases that require modeling of 
relatively large material domains to capture the complete defonnation field, but where 
atomic and subatomic refinement is needed only in very localized regions. Such 
computational issues arise in modeling crack nucleation and propagation, and in 
modeling dislocation formation and interaction. By using coupled models, the size 
limitations of the atomistic simulation can be minimized by embedding an inner region 
where complex dynamic processes and large deformation gradients exist within an outer 
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domain where the deformation gradients are small so that a continuum finite element 
method (FEM) representation of the material becomes appropriate. 

Over the past decade, various methods have been developed to address different 
problems involving atomistically large material domains [1-12]. The most challenging 
problem in developing these coupled methods is the formulation of a seamless 
computational connection along an interface between different material representations. 
A brief review of several representative coupling procedures follows to illustrate the 
current state-of-the-art. 

In coupling atomistic and continuum material representations, the continuity of 
material properties must be maintained while transitioning from individual atoms 
interacting through nonlocal forces to the local stress-strain field fonnalism of continuum 
mechanics. For crack problems, the early efforts of Gumbsch and Beltz [10] led to the 
development of the Finite Element - Atomistic (FEAt) coupling procedure that combined 
an embedded MD system with a finite element domain. A generalized formulation of 
conventional FEM, which allows FEM nodes to be considered as coarse-grained MD 
“atoms” led to another computational scheme for atomistic-continuum coupling called 
Coarse Grained Molecular Dynamics (CGMD). A detailed discussion of CGMD is given 
by Rudd and Broughton [1 1,12]. In yet another coupling method, the Coupling of Length 
Scales (CLS) method [2], the nodes in a finite element model representing the continuum 
region are directly connected to the atoms in an atomistic region fonning an interface of 
“pad” atoms. The region of “pad” atoms, used in this and other atomistic-continuum 
coupling methods, serves to minimize surface tension effects on the atoms in the 
atomistic system but also introduces a constraint due to the elasticity of the interface 
region. The constraining effect of this region is generally considered insignificant and is 
ignored. 

The Quasicontinuum (QC) method, reviewed by Miller and Tadmor [6], is formally 
based on an entirely atomistic description of the material domain. However, for 
computational efficiency, regions are identified in which discrete atoms may be grouped 
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to form a local continuum. The particular representation used is determined by evaluating 
the magnitude of local defonnation gradients and dictates the treatment of 
“representative atoms” or “repatoms.” In the QC formalism, “non-local repatoms” are 
used to represent “real” atoms to form atomistic regions treated by MS/MD methods 
while “local repatoms” are used to define continuum domains by applying both the 
Cauchy-Born rule [13] and aspects of FEM. The interaction of local and nonlocal 
repatoms at the atomistic/continuum interface leads to the generation of “ghost forces” 
that must be mitigated through the introduction of “dead loads” that are iterated for self- 
consistency in the force balance at the interface between subregions. 

Another representative coupling approach is the bridging method of Xiao and 
Belytschko [7] and is based on an overlay approach in which MD and FEM 
representations are superposed in an interface region. This method allows interpolated 
FEM nodal displacements to be associated with atomic displacements in the bridging 
domain. The method explicitly develops coupled energy Hamiltonians for the atomistic 
and continuum regions and enforces compatibility in the bridging domain using Lagrange 
multipliers. This approach also avoids spurious wave reflection at the MD/FEM interface 
without introducing damping or filtering procedures. 

The Coupled Atomistic/Dislocation Dynamics (CADD) method of Shilkrot et al. [8] 
is specifically designed to simulate, identify and pass dislocations between atomistic and 
continuum domains. The method was originally limited to 0 K simulations [8]; but has 
been recently extended to include finite temperature effects in the MD system by linking 
the MD to a quasistatic FEM domain through a thermal damping region [9], Currently, 
CADD uses a two-dimensional material representation due to the complexity of passing 
fully three-dimensional dislocations between the MD and FEM domains. 

A common feature of many of these approaches [1-12] for coupling atomistic and 
continuum representations is the refinement of the finite element (FE) mesh to atomic 
length scales to link the kinematics of the FE nodes to that of the discrete atoms along an 
interface. In this paper, approaches that relate atoms and FE nodes in a one-to-one 
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manner, or through a form of interpolation, will be referred to as direct coupling (DC) 
approaches. 

While DC approaches are straightforward, the fundamental difficulty in their 
development lies in the inherent differences between the atomistic and continuum 
computational models. The physical state of the atomistic region is described through 
nonlocal interatomic forces between discrete atoms of given position and momentum, 
while the physical state of the continuum region is described through continuous stress- 
strain fields that reflect local statistical averages of atomic interactions at larger length 
and time scales. In general, the fonnal connection between continuum and discrete 
quantities can only be achieved through an adequate statistical averaging over scales 
where the discreteness of the atomic structure can be approximated as a continuum. A 
consequence of some DC interfacing strategies in their initial fonnulation, such as QC, 
required that the analysis be performed quasistatic ally at 0 K. Further details of the direct 
coupling methods may be found in the original publications and in several general review 
papers [4-6]. 

In this paper, an alternative to the DC approaches is proposed to construct a coupled 
MD-FEM system. The approach is based on solving a coupled boundary value problem 
(BVP) at the MD/FE interface for a MD region embedded within a FEM domain. The 
method uses statistical averaging over both time and volume in atomistic subdomains at 
the MD/FE interface to determine nodal displacement boundary conditions for the 
continuum FE model. These enforced displacements, in turn, generate interface reaction 
forces that are applied as constant traction boundary conditions [14-16] between updates 
of the FEM solution to the atoms within the localized MD subdomains. Thus, the present 
approach may be described as a local-nonlocal BVP because it relates local continuum 
nodal quantities with nonlocal statistical averages of atomistic quantities over selected 
atomic subdomains. An iterative procedure between the MD statistical displacements and 
the FEM reaction forces ensures continuity at the interface. In this way, the problem of 
redefining continuum variables at the atomic scale is avoided, and the developed 
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interface approach links different time and length scales between the MD region and 
FEM domain. 

With the emphasis of using statistical averages to couple the two computational 
schemes, the developed approach is identified as a statistical coupling (SC) approach. 
Based on the SC approach, the developed MD-FEM coupling method is referred to as the 
embedded statistical coupling method (ESCM). ESCM provides an enhanced coupling 
methodology that is inherently applicable to three-dimensional domains, avoids 
discretization of the continuum model to atomic scale resolution, pennits arbitrary 
temperatures to be applied, and treats, in a rigorous manner, the compensation of surface 
effects in the MD system. 

This paper will detail the ESCM approach for coupling MD and FEM computational 
domains for the case of systems that reach thermodynamic equilibrium or evolve 
quasistatically. While there is no principle difficulty in implementing this approach for 
non-equilibrium systems, it is beneficial to first consider the case of equilibrium 
simulations to illustrate initial applications of this methodology. 

The remainder of this paper is organized as follows. Section 2 describes the structure 
of the coupled MD-FEM model. This includes discussions of the MD and FEM material 
representations, the coupling interface, and the iterative MD-FEM coupling methodology. 
Section 3 presents several validation studies to substantiate the accuracy of the developed 
methodology. Section 4 presents concluding remarks on the overall effectiveness of 
ESCM. Details of internal force calculations involved in the coupling procedure and a 
discussion of model generation are contained in separate appendices. 

2. The ESCM Model 

The ESCM approach is developed to reduce computational costs incurred while 
simulating “large” volumes of material by embedding an inner atomistic MD system 
within a surrounding continuum FEM domain. In principal, the shape of the atomistic 
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region may be arbitrary as shown in Figure 1; however, for simplicity, the special case of 
a circular region is utilized in the present work. Similarly, although any constitutive 
behavior may be assumed for the FEM domain, the present study considers a linear 
elastic continuum. 



5y(x), 5 X y(x) 


>X 


Figure 1. Model topology used in the ESCM approach depicting the continuum 
region represented by the finite element method (FEM), the atomistic domain 
simulated by molecular dynamics (MD), and the MD/FE Interface region coupling 
the two computational approaches. 


The structure of the ESCM model consists of four regions: 1) an Inner MD Region; 2) 
an Interface MD Region wherein MD and finite elements are superimposed; 3) a Surface 
MD Region that does not interact with the FE nodes but is used to compensate for atomic 
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free edge effects; and 4) a FEM domain in which standard finite element equations apply. 
These four regions are depicted in Figure 2. 



Figure 2. Structure of the MD/FE Interface in the ESCM approach. 


Complete details of the ESCM procedure will be presented by discussing general 
aspects of the MD and FEM computational systems, followed by specific details of the 
MD/FE Interface, the Surface MD Region and the MD-FEM coupling procedure. 

2.1 MD and FEM model components 

The Inner MD Region is used to model material phenomena at the atomistic level. 
This Inner MD Region should be large enough to ensure a statistically smooth transition 
from a continuum to an atomistic representation while modeling any of the types of 
processes (e.g., dislocation formation, void nucleation, or crack propagation) that are 



required by the simulation. Together, the Inner, Interface, and Surface MD Regions 
constitute the complete MD system. 

It is important to emphasize that the partitioning of the MD system into different 
regions is not a physical separation of the system. An atom assigned to a particular 
location freely interacts with atoms in its interaction neighborhood that may reside in a 
different region. Thus, the overall simulation is performed using any conventional MD 
technique without any imposition of direct kinematic constraints. The only difference 
between the three MD regions is that, while the atoms in the Inner MD Region are 
subject only to their interatomic forces, the Interface and Surface MD Regions serve the 
added purpose of facilitating the application of external forces involved in the ESCM 
procedure. 

The addition of a FEM domain pennits a large reduction in the computational cost of 
simulations by replacing the atomistic representation with a continuum model in those 
parts of the system where the deformation gradients are small and atomic-level resolution 
is not necessary. The current application uses the FEM domain to simulate an extended 
material model such that the elastic deformation and load transfer due to applied far- Held 
boundary conditions are accurately transferred to the Inner MD Region. The continuum 
field is currently assumed to be static with linear elastic material properties but other 
applications of ESCM might require the incorporation of nonlinear material behavior, 
such as plasticity or general dynamic response, where nonlinear processes generated in 
the Inner MD Region can be propagated into the continuum. 

2.2 MD/FE interface 

The main role of the MD/FE Interface is to provide a computational linkage between 
the MD region and FEM domain. The atoms that surround a given FE node at the 
interface are partitioned to form a cell in the Interface MD Region, called an interface 
volume cell (IVC), as shown in Figure 2. A similar partitioning is also applied to the 
Surface MD Region, forming surface volume cells (SVCs). The IVCs compute averaged 
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MD displacements at their mass center that are then prescribed as displacement boundary 
conditions to the associated interface finite element nodes. The IVCs need not coincide in 
size or shape with the finite element to which the FE node belongs. In the model 
described in this paper, the IVCs are formed through a Voronoi-type construction [17] by 
selecting those atoms with a common closest finite element node. 

Typically, one finite element at the interface encompasses a region of several hundred 
to several thousand atoms. A lower bound for the number of atoms associated with each 
finite element node is determined by the requirement of obtaining a minimally fluctuating 
average of atomic displacements and minimizing the magnitude of generated gradients in 
the MD region bordering the FEM domain. With an effective average at this scale, the 
discreteness of the atomic structure is homogenized enough so that the FEM domain 
responds to the atomistic region as an extension of the continuum. 

During the coupled MD-FEM simulation, a spatial average within each k' h IVC is 
performed to yield the center of mass displacement, 8 ^ k , which is further averaged 

over a certain period of M MD time steps to yield the statistical displacement vector, 

3/,k = (^CM,k ) = T T X (vvw ,/c (tj ) ~ ? C M,k (o)) (1) 

' ' t JVL j = l 

N / \ , 

X 7 t [tj j is the center of mass of the k IVC 

1=1 

containing atoms at positions r. at time tj of the / h MD step. The mass center 

displacement, 8 CMM , in Equation (1) is calculated relative to the initial zero-displacement 
position of the k th IVC, r CM k (0 ) . The determination of this initial position is discussed in 

Appendix A. In turn, the IVCs distribute reaction forces from the interface finite element 
nodes as external forces applied to the corresponding atoms within the IVC. 


In the above expression, 


'CM 


y 


i 
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2.3 Surface MD Region 


In order for the MD domain to deform freely in response to applied reaction forces, it 
is modeled using free surface boundary conditions as discussed in [14,15]. However, the 
existence of a free surface introduces several undesirable effects in the MD system. First, 
it creates surface tension forces that must be removed to avoid distorting the MD 
response. Second, because atoms at or near the free surface do not have a complete set of 
interacting neighboring atoms, the coordination number of the surface atoms is reduced 
so they are less strongly bonded to the surrounding atomic field than those within the 
interior. Under sufficiently large reaction forces, these atoms may be separated from the 
surface layer causing a surface degradation within the MD domain. To mitigate these free 
surface effects and to stabilize the atoms in the Interface MD Region, an additional 
volume of outlying atoms constituting a Surface MD Region is introduced as shown in 
Figure 2. 

While the Surface MD Region eliminates the free surface effects within the Inner MD 
Region, it also introduces an undesirable fictitious stiffness that elastically constrains the 
deformation of the Inner MD Region. The separate effects of surface tension and the 
fictitious stiffness cannot be computed independently. However, their combined effect 

may be defined as a resultant force, f s , which acts at the boundary between the Surface 
MD Region and the Interface MD Region, and is given by the sum of two components 
expressed as 


/,=! + ?• ( 2 ) 

In Equation (2), f is the elastic reaction of the Surface MD Region under deformation, 
and f is the force that results from the surface tension. Both forces are shown in Figure 
3. 

A way to mitigate both the surface tension and the elastic response of the Surface MD 
Region is to estimate and compensate for the force/, . In the ideal case, when f s is fully 
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compensated, the Surface MD Region acts as if it possesses zero stiffness and 
experiences no surface tension, thereby mitigating spurious influences on the Inner MD 
Region. Subdividing the Surface MD Region into a number of SVCs helps to follow the 

variations of f s along the perimeter of the Interface MD Region. For convenience, the 
partitioning of SVCs can be made to follow the IVC partitioning of the Interface MD 
Region. The resultant force is then calculated individually for each SVC. To compensate 

for /' , a counterforce, /'. , is computed along the IVC/SVC interface and then distributed 
over the atoms of each SVC in a similar manner as the nodal reaction forces are applied 
to the IVCs of the Interface MD Region. The calculation of the counterforce, f c , is 
presented in Appendix B. 



Figure 3. Forces within the Interface and Surface MD Regions. R j are 
FEM reaction forces,^ are the resultant surface region forces, and /■ are 
applied compensating forces. 
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2.4 Phonon Damping 


Both the IVCs and SVCs serve the additional role of providing a dissipative damping 
mechanism for phonons propagating into the interface. Potential sources of phonon 
generation are the application of the FEM reaction forces to the IVCs and the resonant 
elastic oscillations in the dynamic MD region. Phonons can also be generated from within 
the Inner MD Region as a result of simulated atomistic processes. In the current 
application of ESCM, these oscillations must be damped in order to achieve equilibrium 
with the static FEM domain. A number of different damping schemes have been 
addressed in the literature. Holian and Ravelo [18], and more recently, Schafer et al. [19], 
found that applying linearly increasing viscous damping to the atoms in a region 
surrounding the center of the MD system can effectively absorb the intense phonon 
waves coming from a propagating crack. In this scheme, a friction force, fj , is given by 

n = -/v (3) 

and is applied to the atoms of the damped region in proportion to the atom’s velocity, v , 
and an appropriately chosen viscous coefficient % [19,20]. The method is efficient and 
simple to implement. Its drawback is that it erroneously decreases the local temperature 
in the damping region resulting in undesirable strain gradients because of thermal 
contraction. 

To avoid disturbing the thennal field, the damping used in the present method is based 
on a modified fonn of Equation (3), where, instead of being proportional to each 
individual atom’s velocity, fj is set proportional to the group velocity of the mass center 
of a certain volume of N atoms. The frictional force, fj , and the group velocity of the 
mass center, v cm , are given by 

1 N 

n = ~Xv cm ; v cm =— Xv- (4) 

/v i = l 
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By controlling the size of this damping region, one can damp phonons of wavelengths 
larger than the volume size while leaving the shorter wavelengths associated with random 
thennal fluctuations unaffected because their contribution to the group velocity averages 
to zero. Phonons introduced by the FE mesh cannot have wavelengths smaller than the 
distance between the interface FE nodes. Thus, it is convenient to choose both the 
Interface and Surface MD Regions as the volumes in which damping is applied. 

To ensure a gradual stepwise increase of the viscous coefficient, %, from zero to % max , 
the Interface and the Surface MD Regions are additionally subdivided into K layers of 
atoms parallel to the interface line. The layer thickness is equal to the maximum 
interaction distance of the applied interatomic potential (0.65 nm). In each layer, % 
increases by a fixed amount A/ starting from A/ for the innermost layer, which neighbors 
the Inner MD region, to KA/ =% ma x for the outermost layer, at the edge of the MD 
domain. The damping volumes used for calculating the group velocity in Equation (4) are 
the cross-sections of the layers with the IVCs and SVCs. 

A discussion of criteria for the selection of the value of the viscous coefficient, %, can 
be found in [19]. In those instances where viscous wave damping is not adequate (e.g., 
collisional dynamics), the more precise non-reflective boundary condition techniques 
discussed in [20-22] may be implemented. 

2.5 MD-FEM coupling 

The MD-FEM coupling in ESCM is achieved through an iterative equilibration 
scheme between the MD region and the FEM domain. In this scheme, iterations begin 
with displacements at the MD/FE Interface that are calculated as statistical averages over 
the atomic positions within each IVC and averaged over the time of the MD analysis. 
These average displacements are then imposed as displacement boundary conditions, 
{ 8 1 } , on the FEM domain. The resulting FEM B VP is then solved to recover new 
interface reaction forces, {Rj}, resulting from the applied interface displacements and 
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any imposed far- Held loading. The new interface reaction forces, { R , }, are then 
distributed to the atoms in the IVCs, thus defining new constant traction boundary 
conditions on the MD system. Between the FEM solution updates, the traction boundary 
conditions are constant and applied to the MD region to ensure that the elastic field from 
the FEM domain is correctly duplicated in the atomistic region. The MD-FEM iteration 
cycle repeats until a stable equilibrium of both displacements and forces between the 
atomistic and continuum material fields is established at the interface. 


The stiffness of the material in the FEM domain is described through a partitioning of 
the global stiffness matrix into a set of stiffness submatrices, \K a p\, with a,/3 = V, F, I 
indicating: V as variables within the interior of the FEM domain; F as far-field variables, 
and / as interface variables. Using these definitions, the static continuum equations of 
state at the n th FEM update at time t n are given by 


\k vv ] [k vf ] 

[Kfv ] I Xff ] 
\k jv ] | K IF ] 


[K n ] 
[KfA 
[K„ ] 


{dy ( { n )} 


i^V ( { n )} 

* l^F ( { n )} 

> = < 

{r f (t n )} 

.[ ('„)}. 




(5) 


To ensure that the FEM domain has the same elastic properties as the MD system, the 
stiffness terms in the [K a p\ submatrices are calculated from the anisotropic elastic 
constants derived from the MD interatomic potential [23] at the same temperature as the 
MD simulation. The FE model is subjected to two types of displacement boundary 
conditions: (1) the far-field displacements {8f}, which define the load over the entire 

coupled MD-FEM system; and (2), the interface displacements {dj } = {sj ,Sj ), 
which represent the deformation response of the MD system at the 1 st , 2 nd , ... , and k th 
I VC. 

The solution for the unknown displacements in the interior of the FEM domain, {8 r }, 
is given by 
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{MO} = !'■ I ({MOH*<4{MO}-PM{S/ <0» 


( 6 ) 


which allows the calculation of the interface reaction forces, {R, (t n )} = ( R ) , Rf R k , ,..) 
of the 1 st , 2 nd , ... , and k th IVC to be obtained from 

{*,('■)}= [ K iv]M n )}+ [^]M„)}+ ( 7a ) 

together with the far- Held forces of constraint 

WO) = W]M,)}+ [*^]{MO}+ [^]{8 / (t„)} (7b) 

The dynamics of an atom i of mass m l> at position r 0> in the embedded MD regions is 
described by Newton’s equations of motion 

nij fj = f ; i e Inner MD Region 

m i r i =f i + Rj / Nj - % v k cm ; i e (lVC\ Interface MD Region (8) 

m i R = fi + fc / N s ~ X vin 1 e ( SVC )k Surface MD Region 

The atoms in the Inner MD Region experience only the atomic force f. = 

j 

resulting from their y' th neighbors. The term f (iJ ' ) is expressed by Equation (B4) as shown 
in Appendix B. The atoms in the interface region, assigned to a given k th IVC, are 
subjected to the additional external force, R k (Equation (7a)), which is distributed over 
the number of contained atoms, N k . The atoms in the Surface MD Region belonging to a 
given k th SVC experience the additional counterforce, ]\ k , which is distributed over the 
N k s atoms contained in their volume. In order to maintain the continuity of forces 
between adjacent cells, the force distribution is interpolated with a linear (or higher order) 
interpolation between atomic positions as a function of each atom’s distance from the 
mass center of its associated cell. In the present study, a simple linear interpolation was 
applied. For the equations governing the Interface and the Surface MD Regions in 
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Equation (8), the viscous friction force, ‘xy k cm > defined in Equation (4), is applied 
uniformly to the atoms contained within the IVCs and SVCs. 

During the MD integration of Equation (8) for a period A?m = MAt, where M is the 
number of time steps and At is the duration of the time step, the new average 
displacements [d I (t n+x )} are computed from Equation (1). The new atomistic 
displacements for the next FEM update at time t n+x =t n +A t M are reapplied in Equations 
(6) and (7a) to calculate the next iterative update of the recovered forces, (f„ +1 )}- 

During the same time interval, the compensation forces {/)(?)} are also evolving through 

Equation (B15) (in Appendix B). The period At M is selected by a determination of the 
convergence rate to a state of dynamic equilibrium between the MD region and the FE 
domain. Applying a suitable damping force (Equation 4) at the MD side of the interface 
ensures a faster convergence rate. The algorithm for the entire coupled simulation is 
summarized in Figure 4. 

Indication for convergence between the MD and FEM domains is the convergence of 
the interface forces {R]} and displacements j 8 [ } to equilibrium steady-state values. This 
convergence can be achieved only if the MD system can reach a dynamic force balance 
with the surrounding FEM system. A convergence criteria will be derived based on the 
force balance between the atomic forces in the MD system and the reaction forces in the 
FEM domain. At equilibrium, the averaged in time motion of the I VC mass centers is 
zero, or 

i^CM ^ = i^CM = 0 ’ ( 9 a ) 

and the change of the averaged total momentum, (Ap k ) , of any k lh IVC due to the FEM 
reaction force for the period of the MD simulation of M time steps is also zero, or 

Ap k = U m i n (C )At = 0 . (9b) 

n=\ i = 1 
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Define coupled MD-FEM Model 


Set initial far-field boundary conditions: 


Set initial interface displacements: 

{8 ,(/„)} 


Solve FEM system and obtain interface reaction forces 


{*,(». )!= [*V Ife- (>, )} +[K, f R if. R (>. )} (7a) 


Perform MD simulation usinq M time steps 


r = t„ +mAt: m = \..M 



Integrate Newton's equation of motion 


m i i'i = fi ; i e Inner MD Region 

m i n =fi+ Rj /N* -x vL >' 1 e ( IVC )k Interface MD Region 

m i n =fi+ fc I N s~% >' 1 e ( Svc )k Surface MD Region 



Calculate new interface displacements 


t n+ \ =t n +MAt 

1 M 

8/(f„ + i) = (s cm)= — Z( F cm(^)-^cm(0)) (!) 
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Performing the same double summation on the second equation in Equation (8) 
results in 


M 


M Nj / N* 


m = 1 z=l 


1 i = 1 



( 10 ) 


which states that, at equilibrium and in accordance with Newton’s 3 ld law, the FEM 
reaction force becomes equal and opposite to the average MD atomic force in the 
corresponding k tb IVC. Equation (10) thus expresses the establishment of static force 
equilibrium between the MD system and the FEM domain and can be used as a 
convergence criteria for the iterative MD-FEM coupling procedure. 


A discussion of practical issues regarding model generation for applying the ESCM 
procedure is presented in Appendix A. 


3. Numerical Verification of ESCM 
3.1 The simulation models 

Four test cases are considered to investigate the behavior of ESCM. First, the 
effectiveness of the application of compensation forces for the mitigation of surface 
tension effects is examined. Second, the dynamic behavior of the MD system is explored 
by varying the rate and sequence of applied external loads. Third, the stress-strain 
continuity between the MD region and the FEM domain is assessed through comparison 
with an exact solution of an elastically deformed plate with a circular hole. Fourth, a 
simulation of the propagation of an edge crack through the FEM domain into the MD 
system is performed to determine the suitability of ESCM for solving problems related to 
crack growth. 

The model geometry used in all of the verification studies is shown in Figure 5. This 
model consists of a circular Inner MD Region of diameter, cImd, that is embedded in a 
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larger exterior square FEM domain of elastic material with a side dimension of cIfe = 
20 cImd- A general discussion of issues related to model construction in applying ESCM is 
presented in Appendix A, while the specifics of the MD and FE models will be discussed 
next. 


T 



h _ 





Figure 5. Model geometry of the test MD-FEM coupled system. 


3.1.1 The MD model 

The material for the simulation models was chosen to be a perfect crystal of 
aluminum. The atomic properties of aluminum were represented by the embedded atom 
method (EAM) potential of Mishin et al. [23], which was fitted to give the correct zero- 
temperature lattice constant, a 0 = 0.405 nm, elastic constants, cohesive energy, vacancy 
formation energy, etc. A 2000 point tabulation of the potential functions in the interval of 
interatomic distances from r mm = 0.25a o (0.1 mn) to the cut-off radius, r c = 1 ,55a 0 (0.628 
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nm) was used. For accurate coupling, material properties of the FEM model are obtained 
from the EAM interatomic potential. 

The interatomic MD potential has a dominant harmonic component that corresponds 
to a local quadratic variation about the equilibrium state which implies a linear variation 
in material elastic properties; far from equilibrium, anharmonic components of the 
potential simulate a nonlinearity in the material elastic response. The nonlinearity is 
difficult to estimate accurately so the FEM calculations used elastic constants that were 
based on the harmonic component only. This choice required the applied strain 
magnitudes to be small (not exceeding 0.5% in this study) so that the anharmonic 
components of the MD potential are small enough to be neglected and the use of linear 
elastic FEM calculations in the MD-FEM coupling procedure is justified. On the other 
hand, using a tabulated potential instead of an analytical expression sets a lower limit on 
the strain values for correctly reproducing the elastic response. When the strain is so 
small that the resulting atomic displacements are smaller than the interval between two 
tabulated points, A r < (r - r mm )/ 2000 , corresponding to a strain of less than 0.05%, then 
the calculated change of the potential energy would depend entirely on the interpolation 
procedure used between the points. To preserve the correct elastic response, the 
interpolation must correctly reproduce the second order derivatives of the potential, 
which means that it has to be based on a second or higher order polynomial function. In 
this study, a second order polynomial interpolation was applied, which, combined with a 
high precision fifth order Gear predictor-corrector integration scheme for calculating the 
atomic trajectories [24], gave errors in estimating the stress in the system under uniaxial 
strain in the range of 0.1 - 0.5% that were less than 0.5 MPa. This is equivalent to 
detecting a change in strain of the order of 0.007%. The high accuracy in strain 
determination allowed accurate measurements of the MD-FEM coupled model at 
relatively small strains of 0.1 to 0.5%. 

The first three test simulations, presented in Sections 3.2.1 to 3.2.3, use a common 
MD system which will be described here. The forth test is performed on a bicrystal MD 
model which will be described in Section 3.2.4. For the first three tests, the MD domain 


21 



is constructed as a circular disk of monocrystalline aluminum with its main 
crystallographic axes [1 0 0], [0 1 0], and [0 0 1] oriented along the x-, y- and z- 
directions, respectively. The MD system is simulated with periodic boundary conditions 
along the z-direction and with free surface boundary conditions along its perimeter. These 
boundary conditions allow the MD domain to deform in an unconstrained manner in the 
x-y plane under the external reaction forces from the FEM domain while maintaining 
constant zero pressure along the z-direction using the Parrinello-Rahman constant- 
pressure simulation technique [25]. Constant temperature is maintained by applying the 
Nose-Hoover thennostat [26] in the Inner MD Region only. 

The thickness of the model along the z-axis is equal to h = 5 a 0 ~ 2.0 nm (Figure 5). 
Though very thin, the MD system mechanically behaves as an infinitely thick plate due to 
the applied periodic boundary conditions along the z-direction. The interference of the 
atoms with their periodic images in the z-direction is prevented by the fact that h is more 
than three times larger than the range of the interatomic potential, r c , thus preserving the 
local three-dimensional characteristics in the system. 

The test simulations were performed at near zero temperature (T= 10 K) to minimize 
the thennal noise, and at room temperature (T = 300 K) to demonstrate ESCM for more 
practically relevant situations in which thermal effects are important. For the models used 
in the present work, effective viscous wave damping in the Interface and the Surface MD 
Regions was achieved by setting % max = 3 eVps/nnr. For this % and a damping layer with 
a width of 0.65 nm, the effective average temperature in the damping volumes decreased 
by only 10% compared to the bulk temperature. 

Initially, the atomistic part of the MD-FEM system was built from a square plate of a 
perfect crystal of aluminum cut along the main crystallographic axes. The dimensions of 
the plate were varied from 22 nm to 164 nm in the x- and y- directions. As described in 
Appendix A, the system was equilibrated under PBC at zero stress and constant 
temperature. After defining the equilibrium zero displacement positions of the mass 
centers f CM (0) of the IVCs of a circular MD domain inside the square plate, the atoms 


22 



outside the prescribed MD domain were removed. A 2 nm thick Surface MD Region was 
defined at the free surface of the MD domain, leaving the rest of the MD domain for the 
inner MD and the interface MD regions. Additional equilibration was performed under 

free surface boundary conditions to estimate the reference forces f s g ={0 , (see Appendix 

B). 

Although the ESCM approach can be applied to an arbitrary interface geometry, the 
circular shape of the MD domain was chosen for several reasons. First, the circular shape 
is the equilibrium shape of a free body with respect to surface tension (when the 
anisotropy of the surface energy is neglected). As already discussed in Section 2, in the 
ESCM approach the MD system represents a free body placed in the elastic field 
generated by the continuum FEM domain. The MD-FEM coupling is more effective if 
the free body is in equilibrium with itself before the coupling is imposed. Such 
equilibrium is not guaranteed for an MD domain of general shape. Starting with a 
nonequilibrium shape of the MD region may cause the applied counterforces to be 
significant compared to the interatomic forces and, thus, may significantly affect the 
atomic state. Second, a circular MD/FE interface is uniform everywhere around the MD 
system, which significantly simplifies the analyses of its behavior and the surface tension 
compensation procedure as described in Section 2.3. Third, the curvature of a circle (or 
sphere) undergoes minimal changes during deformation, compared to any other shape 
because the radius of curvature is constant and equal to the maximum possible value 
everywhere on the surface. Thus, the change of this radius under strain will be minimal. 
Preserving a nearly constant curvature of the free surface during deformation causes the 
surface tension to also remain nearly constant during the simulation thus facilitating the 
application of the surface tension compensation procedure described in Appendix B. 

3.1.2 The FE model 

The elastic continuum region was modeled using 8-node hybrid-stress hexahedral 
finite elements that have a reduced sensitivity to mesh distortion compared to standard 
displacement-based elements, and allow explicit stiffness coefficients to be analytically 
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derived, thereby minimizing their computational requirements [27,28], The elastic 
constants in the material constitutive matrix were derived from the interatomic potential 
for pure aluminum at T = 10 K. The values were averaged for uniaxial stresses from 100 
to 500 MPa, accounting for the nonlinear material properties, as: Cn = 1 12.7 GPa, Cn = 
59.4 Gpa, and C 44 = 30.6 GPa. These values differ by only 3% from the static, zero 
Kelvin elastic constants reported for this potential in [24]. 

The continuum finite element model contains an open inner region of diameter cImd, 
within which the atomistic domain is embedded. Along its perimeter, 80 nodes at z = 
+h/2 and 80 nodes at z = -hi 2 were placed to form 160 FE interface nodes to 
communicate with the embedded MD system. 

The dimensions of the FE mesh, d FE , cImd and h, as shown in Figure 5, were initially 
defined through the proportions of d FE '■ dMD : h = 20 : 1 : 1 . Then, a direct scaling of the 
FE nodal coordinates was performed such that the dimensions <7 md and h matched the 
dimensions of the MD system. Finally, a second scaling of the FEM system was 
performed to preserve the outer dimension ratio d FF : dMD = 20 : 1 . 

3.2 Numerical test results and analyses 

The four test cases selected to interrogate the essential features of ESCM are presented 
in the following sections. Discussions assessing results and details of the analyses are 
included to thoroughly investigate the verification simulations. 

3.2.1 Verification of the surface tension and the Surface MD Region 
stiffness compensation 

The purpose of this simulation is to estimate the magnitude of the surface forces, their 
effect on coupling the two computational domains, and the ability of the compensation 
procedure to mitigate both surface tension effects and the spurious constraint of the 
Surface MD Region stiffness. The simulation is perfonned for the case of a homogeneous 
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perfect crystal of aluminum. Because the effect of the surface tension is expected to be 
relatively weak, the temperature of the simulations was kept at T= 10 K to minimize the 
thennal noise and to increase the sensitivity of the force and pressure calculations. 

The pressure inside the system due to surface tension is defined as the radial 
component of the stress tensor, p s = <j , r , averaged over an isolated MD system with free 
surface boundary conditions in the x- and y-directions, and periodic, zero pressure 
boundary conditions applied in the z-direction. The virial definition of stress [29] inside 
the MD system was used. The surface pressure increases from 10 to 80 MPa with cImd 
decreasing from 164 to 22 nm as presented in Figure 6. The expected dependence of p s on 
surface tension, y (p s = 2yl i I md ) for a cylindrical nanoparticle [30] is well reproduced. 
The surface tension, estimated from the slope of a linear fit to the results in Figure 6 as y 
= 0.9 J/m , is found to be in the range of the calculated surface energies for the 
interatomic potential used (y s = 0.87, 0.943, and 1.006 J/m 2 for (111), (100) and (110) 
surfaces, respectively [23]). While the values for p s (< 0.1 GPa) are relatively small for a 
MD simulation (where typical loads are of the order of 1 GPa), particularly for small MD 
systems, its contribution should not be neglected. 



Figure 6. Dependence of the pressure p s caused by the surface tension of a 
circular MD domain of radius dMD simulated with free surface boundary 
conditions. The slope is proportional to the surface tension y. 
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The effect of applying a counterforce to compensate for the surface tension is shown 
in Figure 7 for the case of duo = 44 nm. When no surface compensation force is applied, 
the internal pressure, p s , gradually increases from zero due to the evolving effect of 
surface tension forces and approaches the value of p s = 41.5 MPa. In contrast, repeating 
the same simulation with the compensation force applied quickly reduces p s to near zero. 
The short initial increase in the value of p s observed during the first 10 ps is a result of 
the iteration procedure (Equation (B15)) for adjusting the counterforce, f c , having not yet 
reached convergence and the compensation being not yet complete. The compensation 
becomes complete approximately 20 ps after the beginning of the simulation. 



Figure 7. Evolution of the internal pressure of a circular MD system 
of dMD = 44 nm with and without surface tension compensation. 


Under deformation, the stiffness of the surface MD region creates an elastic reaction 
force, which adds to the surface tension. As explained in Section 2.3 and in Appendix B, 
the iteration procedure (B 1 5) for f c was constructed to compensate for both effects. The 
efficiency of this compensation for the surface tension and the surface MD region 
stiffness is demonstrated in Figures 7 and 8. The coupled MD-FEM system (Figure 5) is 
subject to uniaxial far-field tensile strain of s = 0.5%. Figure 8 shows the equilibrium 
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radial and tangential displacements of the interface FEM nodes for the MD domain of 
d MD = 44 mn. The nodal displacements with (open symbols) and without (full symbols) 
counterforce compensation are compared to the theoretical profile (solid lines) of a 
homogeneously deformed circle of radius r given as: 


2 


2 


rs. 


= cos' cp-vsin cp 


rs v 


— = (l + v)cos cp sin cp 


( 11 ) 


where v is the Poisson’s ratio, u r is the radial displacement, and u, is the tangential 
displacement. 



Figure 8. Equilibrium radial u r and tangential u { displacements of the interface FEM 
nodes for an MD domain of d MD = 44 mn. The node displacements with (open 
symbols) and without (full symbols) counterforce compensation are compared to the 
theoretical profile (solid lines) of a circle of radius r, subject to a uniform tensile 
strain s. 


When no compensation is applied, the deviations in the displacements compared to 
the theoretical profile due to the surface tension and stiffness of the surface MD region 
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are noticeable, although relatively small. Applying the counterforce reduces these 
deviations to between 1/3 and 1/5 of their uncompensated values. 

Figure 9 shows the combined effect of surface tension and the elastic stiffness of the 
surface MD zone on the equilibrium stress state of the inner MD region for four different 
MD domain sizes ranging from duo = 22 to 164 mn at 0.5% far-field strain. As the 
dynamics of the MD system strongly depends on its size, to make the simulations 
comparable, the time t is rescaled by an estimated relaxation time t 0 for each system size. 
The relaxation time is defined in the standard way as the period required for an 
exponential variable to decrease 1/e (0.386) of its initial value. In Figure 9, the 
exponential variable is the deviation of the current stress <j xx {t ) at time t from its 
equilibrium value o K established at t/t 0 = 5 while the counterforce f (t) was turned off 


• (i2) 


No surface With surface 

compensation compensation 



Figure 9. Evolution of the longitudinal tensile cr xx and transverse a yy stress 
components of the internal stress of the MD system of a coupled MD-FE model 
simulated at 0.5% far field strain. Data for four circular MD systems of diameters 
ranging from 22 to 164 mn is given before (t/t 0 < 5) and after (t/t 0 > 5) surface 
tension compensation (shown separated by the gray band). 
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The relaxation time t 0 was found to increase with increasing size and mass of the MD 
domain. For du d = 22, 44, 84 and 164 mn (containing 47 600, 190 300, 693 400, and 2 
641 000 atoms respectively), a curve fit to Equation (12) gave t 0 = 5.1, 19.5, 61 and 200 
ps, respectively. 

Initially (t/t 0 < 5), no counterforce was applied in the Surface MD Region and a xx 
equilibrated to a lower level compared to the far-field stress, cr G , of a homogeneously 
strained plate. The stress due to the Poisson contraction of the unconstrained boundaries 
(y-±d FE / 2), Oyy, , also deviated from the expected zero level. The reason for these 
deviations is the combined effect of surface tension and the elastic stiffness of the Surface 
MD Region. The deviations of both <j xx and cr vv decrease proportionally to the increase of 
d M D, as expected because of the decreasing surface-to-volume ratio. When the 
counterforce is applied (t/t 0 > 5), the effect of the spurious forces in the Surface MD 
Region for all of the tested MD systems of t/vm from 22 to 164 mn is almost entirely 
negated. 

3.2.2 Simulation of the dynamics of the coupled MD-FEM system 

The dynamic behavior of the coupled MD-FEM system in the simplest case of a 
uniformly loaded homogeneous aluminum plate is depicted in Figure 10. The figure 
presents the stress response of the Inner MD Region to the applied far-field strain. The 
system is the same as that used for the results in Figure 7, with dMD = 44 mn. In this 
numerical test, the prescribed strain of s xx = 0.5% was reached in two ways: first, by an 
instantaneously applied far-field strain of 0.5% at the outer-boundaries of the FEM 
domain, and second, by five consecutive increments of 0.1% each. In both cases, the 
length of the MD iteration simulation was A/m = MAt = 1 ps with M = 500 and At = 2 fs. 

The tensile component of the stress in the MD system, <j xx , converges to nearly the 
same value for both cases shown in Figure 10 and is very close to the far-field stress of 
the FEM domain, <j 0 ( cr xx — > 354 MPa and cr 0 = 358.5 MPa). Similarly, a yy quickly 
relaxes to zero after a temporary jump in response corresponding to each increase in the 
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applied far- Held strain. The test shows that the state of equilibrium, wherein the MD 
system is in force balance with the FEM domain, is obtained with little dependence on 
step size for mono tonic loading up to 0.5% strain. 



Figure 10. Dynamic response of the coupled MD-FEM system under far- field 
homogeneous tensile strain s xx applied either instantaneously (black) or in five 
consecutive increments (gray). 


The source for the small, but systematic deviation (~1 %) of a xx from <j 0 , seen also in 
Figure 9 after the surface compensation has been restored (t/t 0 > 5), is difficult to 
determine as it is almost in the range of the precision limit of the simulation techniques 
used. Most likely, the deviation is due to the anharmonic part of the MD potential, which 
makes the MD domain slightly softer than the linear elastic FEM region. Another cause 
for this deviation is the slightly inexact compensation for the free surface effects. This is 
most likely the case for the smallest system of cImd = 22 nm, shown in Figure 9, where cr xx 
tends to converge below the stress values obtained for the larger systems. The latter effect 
disappears for larger systems where surface effects are weaker and the surface 
compensation is more efficient, and their internal stress converges to the nominal value. 
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After each strain increment of 0.1%, the MD system reached equilibrium with the 
FEM domain within approximately 25 ps, which is consistent with the estimate for t Q = 
19.5 ps in Figure 7 for d md = 44 nm. This time was approximately the same as for the 
instantaneous jump to 0.5%. The relatively large response time observed in Figure 9 
indicates that using a static FE model (Equation (5)) for the continuum part of the system 
is not suitable for simulating processes in the MD system that are evolving faster than t 0 . 

3.2.3 Test for stress-strain continuity over the MD and the FEM regions: 
elastically deformed plate with a circular hole 

To assess the capability of ESCM for generating compatible stress-strain fields in the 
elastic continuum domain and the MD atomistic region, the classic example of a plate 
with a circular hole subjected to uniaxial loading was used. In addition to having an exact 
elasticity solution for the slightly anisotropic material properties used [31,32], this model 
is particularly convenient for two reasons. First, it provides large stress variations (from 
zero to 2.82cr 0 ) around the hole, which can be used to test the continuity of the stress field 
at the MD/FE interface in the case of large stress gradients. Secondly, keeping the peak 
stress, 2.82cr 0 , well below the theoretical elastic limit of the material (recently estimated 
for the applied interatomic potential to be between 3 and 5 GPa [33]) prevents the 
occurrence of any plastic mechanisms in the MD region, which is not addressed in the 
present study, and a static elastic equilibrium can be achieved everywhere in the system. 
The MD domain with g?md = 84 nm used in Section 3.2.1 was implemented as a starting 
model. After equilibration, a circular hole of 40 nm in diameter was introduced in its 
center (see Figure 1 la). 

For comparison, an equivalent fully continuum three-dimensional anisotropic FE 
model was also simulated with a hole of radius 20 nm. This model was generated within 
a square FE mesh of dimension 1.6 pm by 1.6 pm and having the same elastic properties 
(derived from the interatomic potential) as the FE part of the coupled MD-FE model. 
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In both models, the square plate was deformed at 0.5% uniaxial strain along the x- 
direction through displacement-controlled boundary conditions imposed on the outer 
sides of the FEM system 800 nm away from the hole in the MD domain. At this strain, 
the far-field stress, estimated in Figure 9, is a 0 = 358.5 MPa. 

Starting from an undeformed MD region, the coupled MD-FEM iterative simulation 
was performed until equilibrium was established between the MD domain and the FEM 
domain whereby {S, j and \R, } converged to constant values. The stress field for <j xx , (j vv , 

and a xy stress components of the coupled MD-FEM system was calculated and compared 
with the fully continuum FEM solution that was obtained using the ABAQUS software 
package [34]. This comparison is shown in Figure 11. 

The stress profiles for cr xx and <% taken along sections coincident with the x- and y- 
axes passing through the center of the hole, as shown in Figure 11(a), are presented in 
Figures 12(a) and 12(b), respectively. The continuity of the stress profiles is well 
preserved at the MD/FE interface apart from fluctuations that result from chaotic atomic 
thennal vibrations. Additionally, the stress profiles of the coupled model closely follow 
the stress profiles of the fully continuum FEM analysis of the equivalent model. The 
largest discrepancy is seen for the tangential component of the stress ( o\, y along the x- 
profile in Figure 12(a) and a xx along the v-profile in Figure 12(b)) in the MD system, 
which deviates substantially from the continuum prediction when approaching the surface 
of the hole. As shown in Figure 12(b), the continuum solution closely approaches the 
theoretical value of 2.82cr 0 that is calculated for the slightly anisotropic material used. 
One reason for this discrepancy may be the definition of virial stress, which gives poor 
convergence and erroneous estimates near free surfaces [35]. But more likely, it is that 
the continuum model does not correctly account for the nature of the surface tension, 
which results from the occurrence of missing bonds between the atoms at the free 
surface. From the previous analysis (Figure 3), it was found that the normal pressure of a 
free surface with a curvature radius of 40 nm is p s = 45 MPa. This value agrees well with 
the normal component of the stress estimated from the MD simulation ( <j xx along the x- 
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profile in Figure 12(a), and <j yy along the y- pro file in Figure 12(b)) at the edge of the hole, 
where the corresponding FEM solution approaches zero. 
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Figure 11. Simulation of an open hole specimen. Comparison of the stress field for 
cixxj Oyy, and a xy stress components of the coupled MD-FEM system (a, c, e) with a 
full three-dimensional anisotropic continuum FEM solution (b, d, f). 
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Figure 12. Stress profiles in the open hole specimen for a xx and a yy along (a) 
the x-axis and (b) the y-axis scanned through the center of the hole (as given 
in Figure 1 1(a)). The symbols represent the ESCM simulation data, while the 
full lines represent fully continuum FEM results. For reference, the surface 
tension induced pressure p s is also shown. 
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3.3 Example of an edge crack simulation along a grain 
boundary in aluminum 

An important application of the ESCM technique described in this paper is the 
simulation of atomistic processes related to damage. A typical example is the tip of an 
edge crack, where the idealized elastic stress field decreases as l/ Vr and extends to a 
distance r from the tip. The crack tip stress field is much larger than a MD system can 
computationally simulate. A coupled MD-FE model for this example is presented in 
Figure 13, where the MD system is used to atomistically simulate the plastic zone at the 
crack tip, and the FEM domain is used to provide the continuum elastic boundary 
conditions of a far-field tensile strain of s VY = 2% applied along the y-direction. To 
investigate the applicability of ESCM at higher temperatures, the simulation was carried 
out at room temperature (T= 300 K). 

The dimensions of the system are: h = 2.9 nm, cImd = 45 nm, and d\ \_ = 900 nm. The 
MD system represents a bi-crystal with a high-angle grain boundary formed between the 
two crystals along which the edge crack propagates. In the selected coordinate system of 

the model, the orientation of one of the crystals is: (x:[7 7 10]; y:[5 5 7 ]; z: 

[1 1 0]), and the orientation of the other crystal is: (x:[7 7 10]; y:[5 5 7 ]; z: 

[1 1 0 ]). In this way, both crystals are mirror images of each other relative to the 
crystallographic plane {5 5 7 } , which becomes the plane of the grain boundary (GB) 

between them. The GB thus formed is a <1 1 0> ^99 symmetric tilt GB. Crack 
propagation along this GB has been extensively studied by MD simulations [36-38] at a 
cryogenic temperature of T = 100 K. 

In the stresses discussed in [36-38], it was found that, while in one direction the crack 
becomes blunted by deformation twinning, in the opposite direction it propagates in a 
brittle-like manner with very little dislocation emission. The latter direction has been 
chosen as the propagation direction for the edge crack of the simulation in Figure 13. The 
simulation using the ESCM approach perfonned at room temperature showed higher 
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dislocation plasticity than at T = 100 K [36-38]. The remaining problem is how to 
transfer this plasticity to the FEM continuum. Some work related to this issue has been 
started independently by Curtin et al. [3], Shilkrot et al. [8], and Qu et al. [9], where the 
CADD coupling methodology has been developed to follow and transmit dislocations 
between the atomistic and continuum regions. However, that methodology is not 
employed here. 

What is essential in the example shown in Figure 13 is that the ESCM approach 
preserves the continuity of the stress - strain field at the MD-FEM interface even for a 
dynamic problem such as crack propagation simulated at room temperature. The depicted 
profile shows the continuity of stresses across the boundary when the crack speed is slow 
compared to the elastic response time of the system. In the example shown in Figure 13, 
the crack propagation speed was approximately 100 m/s or on the order of 1/30 the speed 
of sound. 


e = 0.02 
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Figure 13. ESCM simulation of an edge crack propagating along a 
<1 1 0>/X99 symmetric tilt grain boundary in Al at 2% far-field uniaxial 
strain showing shear stress a xy . 
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4. Concluding Remarks 


A new statistical approach to couple MD with FEM simulations, denoted the 
embedded statistical coupling method (ESCM), has been developed. The approach is 
based on solving the boundary value problem through an iterative procedure for both MD 
and FEM systems at their common interface. The two systems are simulated 
independently, and they communicate only through their boundary conditions. The FEM 
system is loaded by far- Held loads applied to the external boundaries and along the 
MD/FE interface by nodal displacement boundary conditions that are obtained as 
statistical averages of the atomic positions in the MD system at the mass centers of 
representative interface volume cells (IVCs) associated with each FE node along the 
interface. The MD system, in turn, is simulated under periodically updated constant 
traction boundary conditions that are obtained from the FEM system as reaction forces to 
the MD displacements at the interface. This iterative approach allows the continuity at the 
MD/FE interface to be achieved at different length and time scales inherent to both 
systems without the need of redefining continuum quantities at the discrete atomic scale 
and atomic quantities at the continuum scale. 

Compared to the widely used direct coupling methods, ESCM does not discretize the 
MD/FE interface to atomistic scales but uses a statistical mapping of atomistic behavior 
to a continuum FEM representation. ESCM presents a simple and flexible technique in 
providing elastic boundary conditions for a MD model and eliminates some of the finite- 
size artifacts inherent to a purely atomistic simulation. Because of the statistical 
connection between the MD and the FEM systems, there is no limitation on the 
temperature of the atomistic system as long as the thennal corrections to the elastic 
properties of the continuum system are considered. 

Using static FEM calculations for the continuum part of the system, ESCM shows 
excellent convergence for systems simulated under static elastic equilibrium and 
preserves stress continuity at the MD/FE interface for systems exhibiting relatively slow 
dynamics governed by the MD simulation of the atomistic part of the system. Additional 
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study needs to be performed to determine if the implementation of a dynamic FEM 
simulation can improve the dynamics of the system away from equilibrium. 

In general, the method is relatively easy to implement. Any conventional FEM code, 
including commercial packages, can be used to solve the FEM part of the model. Only 
small modifications to an otherwise conventional MD code are necessary to apply the 
constant tractions to the MD system. 

The verification simulations perfonned in this study demonstrated the effectiveness of 
ESCM to couple atomistic and continuum modeling into a unified multiscale simulation. 
Several idealized test cases were analyzed to interrogate the behavior of ESCM. First, the 
effectiveness of using the Surface MD Region to provide means to emulate infinite 
boundary conditions for the MD system was verified. Second, the dynamic behavior of 
the coupled MD-FEM system was explored demonstrating convergence to the same 
equilibrium state while varying the rate and sequence of applied loads. Third, the stress- 
strain continuity between the MD region and the FEM domain was validated using the 
model of an elastically deformed plate with a circular hole. Finally, simulating the slow 
propagation of an edge crack was performed to demonstrate the overall representational 
capability of the coupled MD-FE model in a system evolving quasi-statically. 
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APPENDIX A 


ESCM Model Construction 

In the ESCM approach, the statistical basis for numerically coupling the different 
computational frameworks provides a much less restrictive set of interfacing 
requirements compared to DC methods and, thus, allows a greater independence in the 
construction of the associated MD and FE models. This aspect of ESCM, however, 
results in additional preparatory work in constructing the coupled model, primarily 
involving the preparation of the initial state of the MD region. A schematic of a MD 
model is depicted in Figure Al. 

Desired MD model to be Superfluous atoms used to apply 

coupled with FEM model. periodic boundary conditions 



Figure Al. Generation of MD model by equilibrating 
within a larger rectangular MD model under PBCs. 


The construction procedure starts with the definition of the shape and size of the MD 
region and FEM domain. The dimensions of the MD system are defined by the 
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dimensions of the combined Inner, Interface, and Surface MD Regions. The dimensions 
of the FEM domain are selected such that the outer boundary defines the desired overall 
material domain and the inner boundary coincides with the IVC mass centers along the 
MD/FE interface. The FE mesh conveniently provides a regular distribution of node 
locations to be used in a Voronoi construction of the IVCs along the MD/FE interface. 
Appropriate interpolation methods (e.g. linear interpolation employing finite element 
shape functions) must be used to accurately map quantities (e.g. interface displacements) 
between the IVC mass centers and the corresponding FE nodes. 

Additionally, it is important in the construction of the ESCM model that the reference 
states of the MD and FEM systems coincide as closely as possible. For a static FEM 
domain, the displacements are zero when the applied forces are zero. For a MD region, 
displacements are statistical quantities that can include some statistical error in their 
estimate that needs to be minimized. 

Preparatory simulations of the MD system involve thermalization, equilibration, and 
the determination of external compensating forces. The calculation of these compensating 
forces is discussed in Appendix B and are required to maintain the initial atomic 
reference states that are necessary for this domain to exhibit the response of an infinite 
system influenced only by external forces when coupled to the FEM computational 
domain. To perfonn the preparatory simulations, an initial MD model of rectangular 
shape is chosen that is large enough such that the geometry of the desired MD region 
(which in the present work is a circular disc) is contained as a subset. This subset can 
subsequently be extracted by removing the atoms outside the desired MD region 
boundaries (see Figure Al). 

To accurately define the zero-displacement reference state of the MD region, the 
system is thermally equilibrated at zero stress and simulated as a constant number of 
atoms, N, constant pressure, P, and constant temperature, T, (NPT ensemble) under 
periodic boundary conditions (PBC) in all three dimensions. PBCs are necessary during 
this step to avoid the presence of a free surface because the surface tension would 
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produce a pressure, p s , on the surface, resulting in erroneous zero-displacement positions 
for the IVC mass centers. 

After achieving equilibrium, an additional MD simulation under PBCs at zero 
pressure and constant temperature is carried out to obtain the statistical time average of 
the mass centers. Larger IVCs and longer initialization times result in smaller statistical 
errors in the reference state because the statistical error of the estimated averages depends 

on the number of atoms, N, per IVC and the time, t, of the simulation as l/^Nt . 
Therefore, the simulation should be carried out long enough such that any systematic 
error is reduced to a level having negligible influence on the coupled simulation. 

Two additional issues must be addressed in the model generation. First, the width of 
the Surface MD Region will generally not fully contain the domain of atoms having 
equilibrium states disrupted due to free surface effects. Therefore, an additional 
simulation will be required to obtain the forces, f s g ={0} , in the reference state as 

explained in Appendix B (Equation (B15)). Second, the elastic anisotropy of the FEM 
domain is a function of the crystallographic orientation of the atomistic microstructure 
and should be derived from the interatomic potential used in the MD system under the 
equivalent thermal and mechanical loading conditions to avoid mismatch of elastic 
properties at the interface. 

The operations discussed in this appendix complete the model construction. This 
model, together with applied far field boundary conditions, is used to start the first MD - 
FEM iteration of the coupled simulation. 

A summary of the individual steps involved in developing the complete ESCM model 
is presented in Figure A2. 
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Define coupled MD-FEM model 


1) Select overall model dimensions. 

2) Select shape and size of the computational partitions: 

a) Define dimension of FEM domain 

b) Define dimension of MD domain 


MD Domain Construction 


Construct MD system of the desired 
microstructure with periodic boundary 
conditions (PBC) containing the MD 
domain as a subset. 


FEM Domain Construction 


Construct FEM mesh that best acco- 
mmodates the shape and microstructure of 
the MD domain. Specify the far-field 
loading conditions. 


MD eauilibration 


Equilibrate the MD system at desired 
temperature and no load under PBC. 


FEM rescalin 


Rescale the FEM mesh according to the 
equilibrated dimensions of the MD 
domain. 


MD - FEM Interface Construction 

1) Superimpose the FEM mesh on the MD domain. 

2) Partition interface MD region. Discretize into interface volume cells (IVCs) and determine 
mass centers to associate with individual FEM nodes. 

3) Partition surface MD region. Discretize into surface volume cells (SYCs). 


MD initialization 

1) MD simulation using PBC’s to estimate the initial zero displacement positions. 

2) Remove all atoms from the MD system exterior to desired MD domain. 


Estimating the reference force state for the Surface MD Region 

(optional) 

1) Make a copy of the MD domain. 

2) Apply a radial force f to the atoms of the Surface MD Region specific for each SVC. (B12) 

3) Equilibrate the copy of the MD domain using f as constant traction boundary conditions. 

4) During equilibration, adjust f so that 8 7 = 0 . 

5) After equilibrium is reached, store f g ={Q} for each SVC and delete the MD copy. (B14) 


Perform coupled MD-FEM simulation 


Figure A2. Flowchart summarizing the MD-FE model construction. 
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APPENDIX B 


Calculation of Compensating Forces 

As discussed in Section 2.3, there are two spurious forces created within the MD 
system in the ESCM approach that have to be eliminated. One is the surface tension 
force, x , created by the applied free surface boundary conditions at the perimeter of the 
MD system. The other is the elastic reaction force, £, , of the Surface MD Region created 
by its stiffness as the MD system deforms. The method to neutralize both of these forces 
is based on applying an external counterforce to the atoms within the Surface MD Region 

/„=-(| + f). (Bl) 

The counterforce is specifically detennined for each SVC and is uniformly distributed 
over the atoms of the SVC so that the total sum of forces in the SVC is zero 

7 c +| + f = 0 (B2) 


Because £ and f cannot be estimated directly during the simulation, the expression 
of the counterforce given in (Equation (Bl)) cannot be used to directly determine f c . 
What can be detennined in the MD simulation, however, is the net spurious force,/,, 
that is generated in the surface MD Region and imposed on the remainder of the MD 
system. As illustrated in Figure 3, the free surface force, f s , for a given SVC is defined 
as the force between the SVC and the Interface MD Region. It is calculated individually 
for each SVC as a sum over the pair forces between atoms (/) of the SVC and all of their 
interaction neighbors (/) lying outside the Surface MD Region 

7,(0= X E/";"(0 (B3) 

ie(SVC) j ^{surf .layer) 
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The pair force between atoms (/) and (j) for a potential based on the embedded atom 
method (EAM) [39] can be expressed as 


7 M (')=- 




dr 


(i.j) 


dr 


(i.j) 


:(i,j) 

JiJ) 


(B4) 


where (f) l '{t) is the potential energy of atom (7) at time t, and r (lJ) =r (,) -r (J) with 

r (i.j) _ | j:(‘,j) | 

To ensure that the equilibrium condition for a perfect crystal lattice is satisfied at zero 
pressure and T = 0 K, the following condition must be met 

(B5) 

A requirement on the division of the Surface MD Region is that each SVC must 
occupy a whole number of lattice unit cells. The reason for this requirement is because 
any resultant force at the atomic scale is a sum of attractive and repulsive forces between 
interacting atoms. The equilibrium condition is satisfied only for the special case of a 
complete, periodic, lattice unit cell. Conversely, equilibrium is not satisfied for individual 
atoms or for arbitrarily defined groups of atoms. 

During simulation with evolving deformations and with the presence of a free 
surface, f s is equal to (after averaging the thermal fluctuations, inherent in each 
atomistically calculated force) the sum of and that part of f , indicated as f s , which is 
contained in the Surface MD Region only 

l = I + (B6) 
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In the above equation, an assumption is made that a very thin Surface MD Region 
of a few nanometers thickness may not contain all the surface tension effects, so that 
the total surface tension force is decomposed into two components 

f = f s +fi (B7) 

where z s is the component that is contained in the Surface MD Region, and z, is the 

component that is contained in the remaining inner part of the MD system. Only f s 
would give contribution to f s in Equation (B6). 

When the counterforce, f c , is applied, /' now equals 

/,=£ + ?, +f c (B8) 

Using the definition for f c expressed by Equation (Bl) for full compensation yields 

fs =I + U “(I + ?)=-?/ (B9) 

Equation (B9) gives the criteria for full compensation of the spurious forces as 
f s = -f , . The counterforce, f c , can now be defined as the force, which is needed to 

maintain f s ( t ) = —z, . This definition has the desirable feature of not requiring that if and 
t be determined explicitly. An iterative procedure is used to maintain f s (t) = -z t by 
correcting fft) within each SVC by the negative of the difference between fft) and 
- Zj at any given time t, as 

X(<+A0V(0-f |7(0+M 7 t (o)= {o} (BIO) 

'‘M 
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Here, At is the MD time step, and /m » At is an adjustable inertial time parameter 
that controls the sensitivity of the counterforce to fast atomic fluctuations of the surface 
force (a suitable choice was found to be t M = 1 OOOAi). 

In order to perfonn the iteration in Equation (BIO), the surface force component, rj, 
has to be determined. If the Surface MD Region is thick enough, f , is a small fraction of 
f , and a good simplifying assumption is to set f , = 0 , which reduces Equation (BIO) to 

X(<+Ar)=7 c (/)-^7,(<) /,(<>)= {0} (Bll) 

Otherwise, when ? 7 is considered non-negligible, the following method can be used 
for its estimation and is based on two considerations. First, when the MD system is not 
deformed, the elastic reaction force of the surface region is zero, or = 0 . Second, the 
assumption is made that the deformation does not change the surface tension, which is 
appropriate if changes in the surface curvature and the surface energy due to defonnation 
are negligible. 

A non-deformed state is defined when the estimated displacements, S , , are zero. 
Here, it is recalled that S, was defined in Equation (1) relative to r CM (0) for a defect free 
system equilibrated under fully 3D periodic boundary conditions with no free surface. To 
make 5 7 = 0 , an external radial force per unit area given by 

f r =-? = -- (B12) 

r 

is unifonnly applied to the atoms of the surface region having a free surface with surface 
tension, y, and radius of curvature, r. Equation (B12) is the Young-Laplace’s equation for 
the internal pressure of a cylindrical particle of radius r due to its surface tension. Even 
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though the Young -Laplace’s equation is strictly applicable to liquids, it has been shown 
that it also provides a reasonable approximation for very small metallic domains [30], 

Under the conditions that 8, = 0 and = 0 , only f s and f r would contribute to f s , 
which can be presented as 

fs = U + fr = t,-t = -?/ (B 13) 

Expressed in another way, recalling Equation (B3), f, is defined as: 

=- 4 , <B14) 

i e (SVC') j g {surf .layer) 

Equation (B14) allows f 7 to be calculated through the atomistic forces in an equilibrated 
MD system when 8 1 = 0 . To avoid uncertainties from random thermal fluctuations in the 
atomistic calculation of f g ={0} , a suitable time averaging at equilibrium can be used in 
Equation (B14). 

After substituting Equation (B14) in Equation (BIO), the iteration procedure for 
calculating the counterforce becomes 

X(<+A«)=/ i (0-^K(0-/,| M o I } £(0)={0} (B15) 

l M 

The iteration in Equation (B15) is performed separately for each SVC to obtain and 
update f] (?) at every MD time step to counteract both f and f at the same time during 
the coupled simulation by maintaining f s (t)~f s g ={0 , . Introducing the initial force state 
f - =!0; as a reference state in Equation (B15) allows the use of a very thin Surface MD 
Region ( 1 to 2 nrn has been used in the present study), so that its stiffness could be small 
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and the compensating force, f c (t), would be a small perturbation to the interatomic 
forces. For a thicker Surface MD Region, r 7 — » 0 , and the simplified iteration (Equation 
(B 1 1)) can be used instead. 
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